Feedback control of unstable steady states of flow past a flat 
plate using reduced-order estimators 

Sunil Ahuja and Clarence W. Rowley 
{sahuja, cwrowley} at princeton.edu 
Mechanical and Aerospace Engineering, 
Princeton University, Princeton, NJ 08544^ USA. 

February 7, 2009 



Abstract 

We present an estimator-based control design procedure for flow control, using 
reduced-order models of the governing equations, linearized about a possibly unstable 
steady state. The reduced-order models are obtained using an approximate balanced 
truncation method that retains the most controllable and observable modes of the sys- 
tem. The original method is valid only for stable linear systems, and in this paper, 
we present an extension to unstable linear systems. The dynamics on the unstable 
subspace are represented by projecting the original equations onto the global unstable 
eigenmodes, assumed to be small in number. A snapshot-based algorithm is devel- 
oped, using approximate balanced truncation, for obtaining a reduced-order model of 
the dynamics on the stable subspace. 

The proposed algorithm is used to study feedback control of two-dimensional flow 
over a flat plate at a low Reynolds number and at large angles of attack, where the 
natural flow is vortex shedding, though there also exists an unstable steady state. For 
control design, we derive reduced-order models valid in the neighborhood of this unstable 
steady state. The actuation is modeled as a localized body force near the leading edge 
of the flat plate, and the sensors are two velocity measurements in the near- wake of 
the plate. A reduced-order Kalman fllter is developed based on these models and is 
shown to accurately reconstruct the flow fleld from the sensor measurements, and the 
resulting estimator-based control is shown to stabilize the unstable steady state. For 
small perturbations of the steady state, the model accurately predicts the response of 
the full simulation. Furthermore, the resulting controller is even able to suppress the 
stable periodic vortex shedding, where the nonlinear effects are strong, thus implying a 
large domain of attraction of the stabilized steady state. 



1 Introduction 



The goal of this paper is two-fold; the first goal is to present an algorithm for develop- 
ing reduced-order models of the input-output dynamics of unstable high-dimensional linear 
state-space systems (such as linearized Navier-Stokes equations with actuation and sens- 
ing), while the second goal is to demonstrate the algorithm by developing estimation-based 
controllers to stabilize unstable steady states of a two dimensional low-Reynolds-number 
flow past a flat plate at a large angle of attack. 



1.1 Model reduction for unstable systems 



Development of feedback control strategies based on linearized Navier-Stokes equations is 
attractive due to the ready availability of a large class of control techniques, and there has 
been substantial progress in this direction in the past decade, reviewed in detail by |Kim fc 



Bewley (2007). However, many of these techniques are limited to relatively small dimen- 
sional systems ~ O(IO^), while the numerical discretization of fluid flows invariably result 
in huge dimensional systems, typically O(10^~^). Thus, model reduction has played an 
important role in making these tools further accessible to fluid flows. 

Extensive research effort in model reduction has focused on the method of proper orthogonal 



decomposition (POD) and Galerkin projection, developed first by Lumley (1970). The main 



disadvantage of this technique is that, although the POD modes capture the energetically 
important structures of the flow, the reduced-order models obtained by the subsequent 
Galerkin projection of the governing equations onto these modes often do not faithfully 
represent the dynamics. Various modifications to improve this method have been proposed 
and used for flow control; refer to the introduction of Siegel et al ([2008) for a review of 
these techniques. The POD/Galerkin methods have been applied for flow control in various 



contexts, such as bluff-body wake suppression (Siegel et al (2008), Graham et al (1999) 



Noack et a/.| ( |2004|), [Tadmor et a/.| ( |2007l )), noise reduction in cavity flow rtRowley &; Jutti- 



judata (2005), Gloerfelt (2008)), and drag reduction in turbulent boundary layers (Lumley 



&: Blossey (1998), Prabhu et al ( |2001[ )). Another model-reduction technique, based on pro- 



jection onto the global (stable or unstable) e igen modes of the flow linear i zed ab out steady 
states, has been used by |Akervik et al (2007) and Henningson k Akervikl tOOSh in the con- 
text of spatially developing flows such as separated boundary layers. In this paper, we focus 
on an approximate balanced truncation method developed by Rowley| ( |2005) as an approx- 
imation to the original method of Moore (1981). This technique captures the dynamically 



important modes of the system, and the non- approximate version provides rigorous bounds 
for the resulting reduced-order models. The method, sometimes called balanced POD, was 



used to obtain models of the linearized channel flow by Ilak &; Rowley (2008) and the Bla- 



sius boundary layer by Bagheri &; Henningson (2009), and was shown to accurately capture 



2 



the control actuation and also to outperform the POD/Galerkin models. 



The balanced truncation method of Moore (1981) is applicable only to systems linearized 



about stable steady states. An extension to unstable linear systems was proposed by Zhou 



et al. (1999), by introducing frequency-domain definitions of controllability and observabil- 



ity Gramians. Reduced-order models were obtained by first decoupling the dynamics on 
the stable and unstable subspaces, and then truncating the relatively uncontrollable and 
unobservable modes on each of the two subspaces. In this paper, we present an approxima- 
tion algorithm for balanced truncation of linear unstable systems, which results in models 
that are equivalent to those of [Zhou et al. (1999) on the stable subspace. The dynamics 
on the unstable subspa ce is treated exactly by a projection onto the global eigenmodes, as 
in lAkervik et al] t007\). 



1.2 Control of flow past 2-D wings 



As a proof-of-concept study, the modeling procedure is applied to the problem of two di- 
mensional low- Reynolds-number flow past a flat plate at a large angle of attack. We develop 
reduced-order models and design controllers that stabilize the unstable steady states of this 
flow. Our motivation for the choice of this problem comes from our interest in regulating 
vortices in separated flows behind low aspect-ratio wings, which is of importance in design 
of micro air vehicles (MAVs). Recently, design of MAVs has been inspired from experimen- 



tal observations in insect and bird flights of a stabilizing leading edge vortex (see Birch &; 



Dickinson (2001) and Ellington et al. (1996j)), which remains attached throughout the wing 
stroke and provides enhanced lift. So, it could be beneficial to design controllers that can 
manipulate the wake of MAVs to enhance lift and achieve better maneuverability in pres- 
ence of wind gusts. Recent studies in this direction, using open-loop control of the flow past 
low-aspect-ratio wings using steady or periodic blowing, were performed computationally 
by 'Taira & Coloniusl (|2009a|) and experimentally by [Williams et al] (|2008|). These studies 



explored different forcing amplitudes and frequencies, locations and directions. However, 
the design of feedback controllers remains a challenge, due to the large dimensionality of the 
problem and complex flow physics. We present computational tools, that we hope can at 
least pave a direction and provide techniques towards addressing some of these challenges. 
We consider the 2-D flow past a flat plate, actuated by a localized body force close to the 
leading edge, with two near-wake velocity sensors. We design a reduced-order compensator 
and show that it is able to suppress vortex shedding at high angles of attack. 

Many previous studies have focused on the control of flow past a cylinder, which is quali- 
tatively similar to the flow past a flat plate at large angle of attack, with the natural flow 
in both the cases being periodic vortex shedding. For the cylinder, the flow undergoes a 
transition from steady state to periodic shedding with increasing Reynolds number, while 
a similar transition occurs in the flat plate with increasing angle of attack. There has been 
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considerable research effort on suppression of this shedding in cyhnder and other bluff body 



wakes, using passive and active, open-loop and feedback control, as reviewed by |Choi et al. 
(2008). Among those, some techniques are based on reduced-order models; for instance. 



Gillies (1998) developed models using artificial neural networks and a POD basis, Graham 



et al] ( |1999| ) modified the POD/Galerkin method to account for actuation by means of 
cylinder-rotation, while Siegel et al] ( |2008[ ) developed a double POD method to account 
for changes in the wake structure during transients. Some earlier efforts in the control of a 
flat-plate wake include those by jCortelezzi (1996), Cortelezzi et al (1997) who used vortex- 
based methods to model the flow past a vertical plate (angle of attack = 90°); vortex-based 



models form their own class of modeling techniques reviewed recently by [Protas (2008). 



Lagrangian coherent structures were used by Wang et al. (2003) to enhance mixing in flow 
past a bluff body with the trailing surface similar to the vertical flat plate. One of the 



few efforts towards control of flat plate at an angle of incidence was by [Zannetti fc lollo 
(2003), who used a passive leading-edge suction control along with a potential flow vortex 



model. Pastoor et al. (2008) also used reduced-order vortex models for drag reduction on 



an elongated D-shaped bluff body. 

This paper is organized as follows: In section[2j we first briefly describe the balanced trunca- 
tion method for unstable systems as developed by Zhou et al. (1999), and the approximate 



balanced truncation procedure called balanced POD of Rowley (2005) for large dimensional 
stable systems. Then, we present an algorithm for approximate balanced truncation of 
large dimensional unstable systems, assuming that the dimension of unstable subspace is 
small and the corresponding global eigenmodes can be computed. In section [3| we briefly 
describe the numerical technique of Colonius &; Taira (2008) using a fast immersed bound- 
ary method, and present the linearized and adjoint formulations of this numerical method. 
In section [4j we present numerical results, using the model example of two-dimensional 
flow past a flat plate at a large angle of attack and a low Reynolds number Re = 100. 
First, we perform a steady- state analysis and compute the branch of steady states in the 
entire range of angles of attack, < a < 90°. We also compute the left and right global 
eigenmodes of the flow linearized about the unstable steady state at o: = 35°. We present 
reduced-order models of the linearized dynamics and use linear optimal control techniques 
to design controllers that stabilize this steady state. Full-state feedback and more practical 
near- wake velocity-measurement based feedback controllers are derived, implemented in the 
nonlinear equations, and shown to suppress vortex shedding. The paper concludes with a 
brief discussion in section [H 
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2 Model reduction methodology 



2.1 Balanced truncation of unstable systems 



We briefly describe a model reduction procedure using the balanced truncation method for 



unstable systems developed by Zhou et al. (1999). Consider the state-space system 



X — Ax + Bu^ 
y = Cx, 



(1) 



where x G is the state, G is the input, and y is the output of the system; the 
dot over x represents differentiation with respect to time. The eigenvalues of A are assumed 
to be anywhere on the complex plane, except on the imaginary axis. 



The standard balanced truncation procedure developed by Moore (1981), valid only for 



stable systems, starts with defining the controllability and observability Gramians of the 
system ([T]) as follows: 







and Wo^ I 'C*Ce'''dt, 
'o 



(2) 



where asterisks are used to denote adjoint operators. A co-ordinate transformation is then 
obtained such that the Gramians ^ of the transformed system are equal and diagonal. 
The diagonal entries of the transformed Gramians, called Hankel singular values (HSVs), 
decrease monotonically and are directly related to the controllability and observability of 
the corresponding states. A reduced-order model is obtained by truncating the states with 
relatively small HSVs, that is, the states which are almost uncontrollable and unobservable. 

For unstable systems, the integrals in ^ are unbounded and hence the Gramians are ill- 
defined. A modified technique was proposed by Zhou et al (1999) based on the following 
frequency-domain definitions of the Gramians: 



- / {jul - A)-^BB%-jcjI - dcj, 

J-oo 
o 

{-jul - A*)-^C*C{juI - A)-^ duj. 



27r 
1 

2^ 



(3) 
(4) 



By using Parseval's theorem, it can be shown that for stable systems, the frequency- domain 
definitions ^ are equiv alent to the time-domain definitions The model-reduction 

procedure of Zhou et al. (1999) begins by first transforming the system ([T]) to coordinates 



5 



in which the stable and unstable dynamics are decoupled. That is, let T be a transformation 
such that if X = Tx, the system ([T]) transforms to 

y = {Cu Cs) X. (5) 

Here, Au and As are such that all their eigenvalues are in the right- and left-half complex 
planes respectively, while Xu and Xg are the corresponding states. Next, denote the con- 
trollability and observability Gramians corresponding to the set (Ag, Bg^Cg) describing the 
stable dynamics by and respectively. Similarly, denote the Gramians corresponding 
to the set {—Au,Bu,Cu) by and W^. The Gramians of the original system are then 
related to those corresponding to the subsystems by: 

W^J 

and Wo = {T-'r(^^^ Wll^~'- 

A system is said to be balanced if its Gramians defined by ^ are equal and diagonal, 
in which case the diagonal entries are called the generalized Hankel singular values. A 
reduced-order model is obtained by truncating the states with small generalized HSVs. 



2.2 Approximate balanced truncation of stable systems 

For systems of large dimension ^ O(10^~^) such as those encountered here, the Grami- 
ans ([g]) are huge matrices which cannot be easily computed or stored. A computationally 



tractable procedure was introduced by Rowley (2005) for obtaining an approximate bal- 



ancing transformation. We first briefly describe this method, valid only for stable systems, 
and then present an extension to unstable systems. The procedure consists of computing 
the impulse responses of the system ([T]) and stacking the resulting snapshots of the state x 
as columns of a matrix X. It also requires state-snapshots of the impulse responses of the 
adjoint system 

w = (7) 

which are stacked as columns of a matrix Z. Then, the Gramians ^ can be approximated 
as 

Wc^XX"", Wo^ZZ\ (8) 
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The approximate Gramians ([8]) are not actually computed due to the large storage cost, 
but the leading columns (or modes) of the transformation that balances these Gramians 
are computed using a cost-efficient algorithm. It involves computing the singular value 
decomposition of 

Z-X = U^V'=(U^ U,)(^^ (9) 

where Si G M^^^ is a diagonal matrix of the most significant HSVs greater than a cut-off 
which is a modeling parameter, while XI2 G is a diagonal matrix of smaller and 

zero HSVs. Note that Z*X G R^«><^^ is a small matrix, where rig and rio are the number 
of snapshots of the impulse responses of systems ([T]) and Q respectively. For fluid systems 
that we are interested in, the typical number of snapshots is O(10^~^), thus resulting in a 
reasonable computational cost, and typically r < 50. The leading columns and rows of the 
balancing transformation and its inverse are obtained using: 

$ = XViY.~^'^, * = ZUiT.~^'^, (10) 

where G R^^^, and the two sets of modes are bi-orthogonal; that is, ^^*$ = /. The 
reduced-order model of ([T]) is then obtained by expressing x — $a, a G M^, and using the 
bi-orthogonality of $ and 

d = ^*A$a + **Bi^, (11) 
y = C$a. (12) 



2.2.1 Output projection 

When the number of outputs of the system (rows of C) is large, the algorithm described in 
section 12.21 can become intractable. The reason for this is that it involves one simulation 
of the adjoint system Q for each component of the dimension of which is the same 
as the number of outputs. This number is often large in fluids systems where a good 
model needs to capture the response of the entire system {C = /) to a given input. To 



overcome this problem, Rowley (2005) proposed a technique called output projection^ which 
involves projecting the output y of ([T]) onto a small number of energetically important modes 
obtained using proper orthogonal decomposition (POD). Let the columns of O G M^^^ 
consist of the first m POD modes of the dataset consisting of outputs obtained from an 
impulse response of ([T]). Then, for the purpose of obtaining a reduced-order model, the 
output of ([1]) is approximated by 

y = ee*Cx, (13) 

where 00* is an orthogonal projection of the output onto the first m POD modes. The 
resulting output-projected system is optimally close (in the L^— sense) to the original sys- 
tem, for an output of fixed rank m. With this approximation, only m adjoint simulations 
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are required to approximate the observability Gramian; refer to Rowley (2005) for details. 
The number of POD modes m for output projection is a design parameter and is typically 
chosen to capture more than 90% of the output energy. In the rest of this paper, the models 
resulting from this approximation of the output are referred to as m—mode output projected 
models. 



2.3 Approximate balanced truncation of unstable systems 



The approximate balancing procedure described in the previous section, which is essentially 
a snapshot-based method, does not extend to unstable systems since the impulse responses 
of Q and are unbounded. We could consider applying the algorithm to the two sub- 
systems given in ([5]), but the transformation T that decouples ([T]) itself is not available. 
However, when the dimension of the unstable sub-system is small, we show that it is not 
necessary to compute the entire transformation T and it is still possible to obtain an ap- 
proximate balancing transformation. Here, we present an algorithm for computing such a 
transformation and also show that it essentially results in a method that is a slight variant of 
the technique of Zhou et al. ( 1999[ ) presented in section 2.1, The idea behind the algorithm 
is to project the original system ([T]) onto the stable subspace of A. Then, one obtains a 
reduced-order model of the projected system using the snapshot-based procedure described 
in section [2^ The dynamics projected onto the unstable subspace can be treated exactly 



on account of its low dimensionality. 

We assume that the number of unstable eigenvalues Uu is 0(10) and can be computed numer- 



ically, say using the computational package ARPACK developed by Lehoucq et al (1998) 



We further assume that the bases for the right and the left unstable eigenspaces ^u^'^u ^ 
^nxnu f^^^ computed. For the algorithm, we need the following projection operator onto 
the stable subspace of A: 



(14) 



where and have been scaled such that = luu- We use the operator Pg to obtain 

the dynamics of ([T]) restricted to the stable subspace of A as follows: 



Xs = Axs + f'sBu, 

Us — ClPs^s 



(15) 



where The adjoint of ([l5| is the same as the dynamics of ^ restricted to the 

stable subspace of A* using P* , and is given by 



A*Zs + ¥lC*v, 



(16) 
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where Zs = P*^:. We compute the state-impulse responses of (15) and JTg]) and stack the 
resulting snapshots Xs and Zs in matrices Xs and Zs respectively. As in ([9]), we compute the 
singular valued decomposition of Z^Xs and use the expressions (10) to obtain the balancing 



modes $5 and the adjoint modes ^ si where again = Ir- The reduced-order modes are 

obtained by expressing the state X as 



X 



(17) 



where au G 
we obtain 



and as . Substituting ([T7| in ([T]) and pre- multiplying by and ^^*, 



da 



d 
dt 



y 



^ u 



+ 



Bu 



a. 



(18) 
(19) 



Now, since TQjige{A^u) ^ span($^i), we can write A^^ 
using the properties of eigenvectors, we have ^"lA^u - 
shown that ^lA^ 



= $^iA for some A G R^^^^^, and 
^^*$^iA = 0. Similarly, it can be 
0. Thus, the cross terms in (18) are zero and the reduced-order model 



IS 



da ^ f^lA^u 
dt \ 



+ 



Bu 




(20) 



The procedure described so far to obtain the reduced-order model (20) is related to the pro- 



cedure of [Zhou et a/.| ( |1999| ) described in section |2.1[ It can be shown that the transformation 
that balances the Gramians defined by ^ results in a system in which the unstable and 
stable dynamics are decoupled. Furthermore, the resulting stable dynamics are the same 
as those given by the equations describing the a^-dynamics of ( [20| ). That is, balancing the 
stable part of the Gramians Wc and Wq defined in ^ (balancing and W^) is the same 
as balancing the Gramians of the stable subsystem ( |l5[ ); a proof is outlined in appendix [A} 
In our algorithm, the unstable dynamics are not balanced, while they are by Zhou et al. 



(1999). A disadvantage of Zhou's approach is that an unstable mode might be truncated 



resulting in a model which does not capture the instability, which is undesirable for control 
purposes. 



2.3.1 Output projection for the stable subspace 



For systems with a large number of outputs, the number of adjoint simulations (16) can 
become intractable; however, the output projection of section [2.2. 1| can readily be extended 
to unstable systems. Instead of projecting the entire output y onto POD modes, we first 
express the state x — Xu + Xs^ where Xu — {I — '^s)x and Xs — P5X are projections on the 
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unstable and stable subspaces of A respectively. We similarly express the output as y = 
Uu -\- Us — C{xu + Xs). We then project the component i/s onto a small number of POD 
modes, of the data set consisting of the outputs from an impulse response of ( p!5| ). If the 
POD modes are represented as columns of the matrix Qg ^ M^^^, the output of ([l]) is 
approximated by 



y = [Cil - Fs 



X 



(21) 



Finally, with the state x expressed by the modal expansion (17), the output of the reduced- 



order model (20) is given by 



(22) 



2.3.2 Algorithm 

The steps involved in obtaining reduced-order models of ([!]), for the case where the output 
is the entire state (C = /), can now be summarized as follows: 

1. Project the original system ([T]) onto the subspace spanned by the stable eigenvectors 
of A in the direction of the unstable eigenvectors of A to obtain (15). Compute the 



(state) response to an impulse on each input of (15) and stack the snapshots in a 
matrix Xg. 

2. Assemble the resulting snapshots, and compute the POD modes 0j of the resulting 
data-set. These POD modes are stacked as columns of O^. 



3. Choose the number of POD modes one wants to use to describe the output of (15). 
For instance, if 10% error is acceptable, and the first m POD modes capture 90% 
of the energy, then the output is the velocity field projected onto the first m modes. 
Thus, the output is represented as i/s = 0*Xs. 

4. Project the adjoint system ([7| onto the subspace spanned by the stable eigenvectors 
of A* in the direction of the unstable eigenvectors of A* to obtain (16). Compute the 



(state) response of (16), starting with each POD mode 9j as the initial condition (one 



simulation for each of the first m modes). Stack the snapshots in a matrix Zg. 

5. Compute the singular value decomposition of M = Z*A* = UsT^sVs^ where G 
R^><^ and r = rank(M). 



6. Define balancing modes ipj and the corresponding adjoint modes tpj as columns of the 

^here 

= XsVs^i^/^ = YsUsT.-^'^. (23) 



matrices $s and ^>a, where 
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7. Obtain the reduced-order model using ( 20 ) , which can be written as 



da 
dt 

y 



Au 
A., 



a + 



Bu 
B.s 



Cu Cg J = Cci 



u = Aa + Bu^ 
where, 



Au = K^^u, Bu = KB, Cu 



B. 



KB, Cs = e.e:$.. 



(24) 

(25) 
(26) 

(27) 
(28) 



Alternatively, the outputs can be considered to be simply the coefficients of the un- 
stable modes and the coefficients of the POD modes 6^ of the stable subspace. 
With this choice, the output can be represented as 



y = 



Cu 
, Cs, 



Ca, where, 



Cu^in^, Cs^ei^s. 



(29) 
(30) 



Finally, if the initial state xq is known, the initial condition of (24) can be obtained 
using 



(31) 



In the remainder of this paper, we demonstrate the algorithm developed in this section 
by obtaining reduced order models of the 2-D uniform flow past a flat plate, and develop 
controllers based on these models to stabilize the unstable steady states that exist at high 
angles of attack. 



3 Numerical scheme 



The numerical scheme used is a fast immersed boundary method developed by [Colonius fc 



Taira (2008), and is briefly described here. The method is then used to develop the linearized 



and adjoint formulations. Consider the following form of the incompressible Navier- Stokes 
equations, based on the continuous analog of the immersed boundary formulation introduced 
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by|Peskiii|([l972D: 



<0 



dtu + u-Vu = -Vp + -^"^^u + J - x)d^. (32) 

V • = 0, (33) 

J u{x^t)5{x — £^)dx — ub^ (34) 

where p and / are the appropriately non-dimensionahzed fluid velocity, pressure and 
surface force respectively. The force / acts as a Lagrange multiplier that imposes the no- 
slip boundary condition on the Lagrangian points ^, which arise from the discretization of 
a body moving with velocity ub- We consider the body to be a stationary flat plate at an 
angle of attack a; that is, here ub = 0- The Reynolds number is deflned as Re = Uc/u 
where U is the free-stream velocity, c is the chord-length and u is the kinematic viscosity. 
Equations ([32j[34]) are discretized in space using a second-order flnite-volume scheme on a 



staggered grid, and a discrete curl operation C^{-) =^ V x (•) is introduced to eliminate the 
pressure and obtain a semi-discrete formulation in terms of the circulation 7: 



^ + C^E^f = -f3C''C-f + C''Af{q) + 6c^, (35) 
ECs = UB = 0. (36) 



The incompressibility condition (33) is implicitly satisfled by an appropriate construction 
of C. The discrete Laplacian is represented by —C^Cj^ using the identity V^7 = V(V • 
7) — V X (V X 7) = —V X (V X 7); the constant f3 = where A is the uniform grid 

spacing. The operator E^ smears the Dirac delta function of (32) over a few grid points. 



The nonlinear term Af{q) is the spatial discretization of g x 7, where q is the discrete velocity 
flux, in turn related to the discrete stream function s and circulation 7 as: 

q^Cs, 7 = C^q, and s = (C^C)-^. (37) 
A uniform grid and a choice of simple boundary conditions result in a fast algorithm. 



With a uniform grid, the discrete Poisson equation (37) is solved by means of the eflicient 
discrete sine transform. The boundary conditions specifled are Dirichlet and Neumann for 
the velocity components normal and tangential to the domain boundaries, which for the 
flow past a flat plate imply a uniform-flow in the far-fleld. These boundary conditions 
are however valid for only a sufliciently large domain, and are imposed by employing a 
computationally eflicient multi-domain approach. The domain is considered to be embedded 
in a series of domains, each twice-as-large as the preceding, with a uniform but a coarser 
grid having the same number of grid points. The Poisson equation, with zero boundary 
conditions, is solved on the largest domain and the stream function is interpolated on the 
boundary of the smaller domain, which are in turn used to solve the Poisson equation on 
the smaller domain. For the flow past a flat plate considered here, the typical size of the 
largest domain is around 40 chord lengths in each direction. Finally, the time-integration is 
performed using the implicit Crank-Nicolson scheme for the viscous terms and the second- 
order accurate, Adams-Bashforth scheme for the convective terms. 
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3.1 Linearized and adjoint equations 



For deriving reduced-order models useful for control design, we first linearize equation (35) 
about a pre-computed steady state (70, qo). The linearized equations are the same as 
equations ( 35|36 ) with the nonlinear term J\f{q) replaced by its linearization about the 



steady state, and is denoted by A/l(7o)7 = ^0 x 7 + ^ x 7o where the flux q is related to 7 



by (37). Thus, the linearized equations are: 

+ CTE^f = -f3C^C^ + C^AAl(7o)7, (38) 



dt 

ECs = 0. (39) 

The boundary conditions for the linearized equations are hc^ = on the outer boundary of 
the largest computational domain. 



In order to derive the reduced-order models using the procedure described earlier, we need 
to perform adjoint simulations. In order to derive the adjoint equations, we define the 
following inner-product on the state-space: 



(71,72) 



/ 7i(C^C)-S2dx. (40) 



That is, the inner-product defined on the state-space is the standard L^-inner product 
weighted with the inverse-Laplacian operator. This choice is convenient as it results in the 
adjoint equations which differ from the linearized equations only in the nonlinear term. A 
derivation is outlined in appendix [B] and the resulting equations are: 

^ + C^E^^ = -PC^CC + {C^C)ML{l^fqa. (41) 
ECi = 0, (42) 

where the variables ^ and ifj are the duals of the discrete circulation 7, stream function s 
and body force /, respectively, and qa — is the dual of flux q. The adjoint of the linearized 
nonlinear term is {C^ C)J\fL{l{))^ qa^ which can be shown to be a spatial discretization 



of V X (70 X qa) — V^(go X qa)- Sincc equation (41) differs from (38) only in the last term 
on the right hand side, the numerical integrator for the adjoint equations can be obtained 
by a small modification to the linearized equations solver. 



The nature of the multi-domain scheme used to approximate the boundary conditions of 
the smallest computational domain, results in a multi-domain discrete Laplacian which is 
not exactly self-adjoint to numerical precision. As a result, the adjoint formulation given 
which also uses the same multi-domain approach, is not precise and results in 
small, rather insignificant, errors in the computation of the reduced-order models. 
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Figure 1: Actuator modeled as a localized body force near the leading edge of the flat plate, 
with the angle of attack fixed at a = 35°. Vorticity contours are plotted, with negative 
contours shown by dashed lines. The velocity-sensor locations are marked by solid circles. 



4 Results: flow past a flat plate 



We apply the model reduction techniques developed in the previous sections to the uniform 
flow past a flat plate in two spatial dimensions, at a low Reynolds number. Re = 100. We 
obtain reduced-order models of a system actuated by means of a localized body force near 
the leading edge of the flat plate; the vorticity contours of the flow field obtained on an 
impulsive input to the actuator are shown in Fig. [TJ Using these reduced-order models, we 
develop feedback controllers that stabilize the unstable steady state at high angles of attack. 
We first assume full-state feedback, but use output projection described in section [2.2.1| to 
considerably decrease the number of outputs in order to make the model computation 
tractable. Later, we relax the full-state feedback assumption, and develop a more practical 
observer-based controller which uses a few velocity measurements in the near-wake of the 
flat plate (shown in Fig. [T]) to reconstruct the entire flow. 



4.1 Numerical parameters 

The grid size used is 250 x 250, with the smallest computational domain given by [—2, 3] x 
[—2.5,2.5], where lengths are non-dimensionalized by the chord of the flat plate, with its 
center located at the origin. We use 5 domains in the multiple-grid scheme, resulting in 
an effective computational domain 2^ times larger the size of the smallest domain; thus 
the largest domain is given by [—32,48] x [—40,40]. The timestep used for all simulations 
is St = 0.01. 
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4.2 Steady-state analysis 



Since our approach is to obtain reduced-order models of the flow linearized about a given 
steady state, we first need to compute these steady states. The model-reduction of unstable 
systems involves projecting the dynamics onto a stable subspace, for which we also need to 
compute the right and left eigenvectors of the linearized dynamics. This section concerns 



this steady-state analysis, using a "timestepper-based" approach as outlined in Tuckerman 



& Barkley (1999 


) and 


Kelley et al. ( 


2004) 



A simple way of computing stable steady states is by simply evolving the time-accurate sim- 
ulation to stationarity. However, unstable steady states cannot be found in this manner, and 
stable steady states near a bifurcation point could take very long to converge. Instead, we 
use a timestepper-based approach which involves writing a computational wrapper around 
the original computational routine to compute the steady states using a Newton iteration. 
If the numerical timestepper advances a circulation field 7^ at a timestep A: to a circulation 
field 7^+^ = $t(7^) after T timesteps, the steady state is given by the field 70 that satisfies 



^(70) = 70 - ^t(7o) 0. 



(43) 



The steady states are given by zeros of ^(70), which could, in principle, be solved for 
using Newton's method. However, the standard Newton's method involves computing and 
inverting Jacobian matrices at each iteration, which is computationally infeasible due to the 
large dimension of fluid systems. Instead of computing the Jacobian, we use a Krylov-space 
based iterative solver called Generalized Minimal Residual Method (GMRES) developed 
by Saad & Schultz (1986) to compute the Newton update (see Kelley ( |1995[ ) and [Trefethen 
& Bau (1997) for a description of the method). This method requires computation of 



only Jacobian- vector products Dg{'y) -v^ which can be approximated using finite differences 
as [^(7+^1;)— ^(7)]/e, for < e ^ 1. So, the Jacobian-vector products can also be computed 
by invoking the appropriately-initialized timestepper. A nice feature of GMRES is relatively 
fast convergence to the steady state when the eigenvalues of the Jacobian Dg('yo) occur in 



clusters; see Kelley (1995) and Kelley et al (2004) for details. For systems with multiple 



time-scales, such as Navier-Stokes, most of the eigenvalues of the continuous Jacobian lie in 
the far-left-half of the complex plane. Thus, the corresponding eigenvalues of the discrete 
Jacobian D^Ti ^ot a sufficiently large value of T, cluster near the origin. 

The procedure described above is used to compute the branch of steady states for the angles 
of attack < a < 90°; the parameter T in ([43]) is fixed to 50 timesteps. The lift and drag 
coefficients, Cl and Cd, and their ratio Cl/Cd with changing a are plotted in F ig, [g} As 
with flow past bluff bodies with increasing Reynolds number (for example, see Provansal 



et al. (1987)), the flow undergoes a Hopf bifurcation from a steady flow to periodic vortex 
shedding as the angle of attack a is increased beyond a critical value ac, which in our 
computations is ac ^ 27°. Also plotted in the figure are the maximum, minimum, and 
mean values of the forces during shedding for a > ac- We see that the (unstable) steady 
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Figure 2: Forces on a flat plate at a fixed angle of attack a and at i?e = 100, showing a 
transition from a stable equilibrium to periodic vortex shedding at a ^ 26. Shown are the 

force coefficients corresponding to the stable ( ) and unstable ( ) steady states, and 

the maximum and minimum ( ), and the mean ( ) values during periodic vortex 

shedding. Also shown are the vorticity contours (negative values in dashed lines) of steady 
states at a = 15°, 55° and the flow fields corresponding to the maximum and minimum 
force coefficients at o: = 55°. 
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Figure 4: Basis vectors of the unstable eigenspace of the hnearized (left) and the adjoint 
(right) equations. Vorticity contours are plotted (negative contours are dashed). 



state values of the lift coefficient are smaller than the minimum for the periodic shedding 
till a ^ 75°, after which they are slightly higher, but still smaller than the mean lift for 
the periodic shedding. The (unstable) steady state drag is much lower than the minimum 
value for periodic shedding. The ratio Cl/Cd of the (unstable) steady state is close to the 
mean value for shedding. Thus, if the large fluctuations in the forces are undesirable at 
high angles of attack, it would be useful to stabilize the unstable state. The steady state 
at q: = 35° is shown in Fig. [sj^a), and a time history of the lift coefficient Cl with this 
steady state as an initial condition is shown in Fig. |3]^b). Since the steady state is unstable, 
the numerical perturbations excite the instability, and the flow eventually transitions to 
periodic vortex shedding. 

We also compute a basis spanning the right and left unstable eigenspaces {^u and ^u) of the 
flow linearized about the unstable steady states, which are required in our model reduction 
procedure, for restricting dynamics onto the stable subspace. As the flow undergoes a Hopf 
bifurcation, a complex pair of eigenvalues crosses the imaginary axis from the left half of 
the complex plane; thus the dimension of the unstable subspace is two. Here, we obtain 
the basis spanning this subspace by a numerical implementation of the power method (see 
page 191 of Trefethen &; Bau ( 1997[ )). We begin the linearized simulation with a very small 



random noise to excite the instability. After a sufficiently long time, the stable modes decay, 
and the dynamics lies close to the unstable subspace. Any two independent snapshots of 
the remaining dynamics gives a good approximation of the basis spanning the unstable 
eigenspace. Similarly, a basis spanning the left unstable eigenspace can be computed by 
initializing the adjoint equations with a small random noise. A basis spanning the right 
and left unstable eigenspaces of the flow linearized about the steady state at a = 35° is 
plotted in Fig. |4| These modes are qualitatively similar to the structures during periodic 
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Figure 5: Vorticity contours of (a) the flow field shown in Fig. [T| projected onto the stable 
subspace, and (b,c) the first- and fifth-most energetic POD modes of the impulse response, 
restricted to the stable subspace. 



vortex shedding, but have different spatial wavelengths, as reported in earlier studies by 



Noack et al (2003) and Barkley (2006) 



4.3 Reduced-order models 



We now describe the process involved in deriving reduced-order models of the input-output 
response ofQ, which in this example are the linearized incompressible Navier-Stokes equa- 



tions ([38| [39|). The actuator used is a locahzed body force close to the leading edge of 
the flat plate, plotted in Fig. [T] The models are derived using the procedure outlined in 
As seen in equation 



section 



2.3 



( |20[ ), the output of the system is considered to be the entire 
velocity field, observed as a projection onto (a) the unstable eigenspace, and (b) the span 
of the leading POD modes of the impulse response restricted to the stable subspace. 

The first step in computing the reduced-order models is to project the flow field B onto 



the stable subspace of (38, 39) using the projection operator defined in equation (14) 



the unstable eigenvectors computed in section |4.2| are used to define Pg numerically. The 
vorticity contours of the corresponding flow field FgB are plotted in Fig.lsk. The next step 



is to compute the impulse response of (15). Instead, for practical reasons, we compute the 
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Figure 6: Vorticity contours of the leading (in the order of Hankel singular values of the 
stable subspace dynamics) first and third balancing (left) and adjoint (right) modes. 



impulse response of 



(44) 



that is, at each timestep of integration, we project the state Xg onto the stable subspace 



of A. Because the stable subspace is an invariant subspace for the linearized dynamics (38) 



theoretically, the impulse responses of equations (15) and (44) are exactly the same, and 



they are the same as that obtained by restricting the impulse response of ([!]) to its stable 
subspace. However, due to the (small) numerical inaccuracy of the projection (which is 
a result of the numerical inaccuracy of the unstable eigenspaces and 'ifu)^ the dynamics 



of (15) is not strictly restricted to the stable subspace and, in the long term, grows without 
bound in the unstable direction. Next, we compute the POD modes 91 of the impulse 



response of (44), and consider the output of (44) to be the state Xg projected onto a certain 



number of these POD modes. Here, 200 snapshots spaced every 50 timesteps were used to 
compute the POD modes. The leading 4 and 10 POD modes contain 84.47% and 99.98% 



of the energy respectively and, as it has been observed in previous studies (see Deane et al 



(1991) and |Ilak fc Rowley (2008)), these modes come in pairs in terms of their energy 
content, a characteristic of traveling structures; the leading first and third POD modes are 
shown in Fig. [5| 

The next step is to compute the adjoint snapshots, with the POD modes of the impulse 
response (projected onto the stable subspace of the adjoint) as the initial conditions. As 
the linearized impulse response, these simulations are also restricted to the stable subspace. 



Again, instead of computing the response of (16), we compute that of the following system: 



(45) 
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Figure 7: The empirical Hankel singular values (- 
controllability ( , o) and observability ( , >< 



and the diagonal elements of the 



-, x) Gramians of a 25-mode model with a 
4, 10, and 20-mode output projection, for the unstable steady state at a = 35. 



The snapshots of the impulse responses of systems ( |44[ ) and ( |45| ) are stacked as columns 
of X and Z, and using the expressions g and we obtain the balancing modes 0^ and 
the adjoint modes tpl- We used 200 snapshots of the linearized simulation and 200 snapshots 
of each adjoint simulation, with the spacing between snapshots fixed to 50 timesteps, to 
compute the balancing transformation. These number of snapshots and the spacing were 
sufficient to accurately compute the modes; further reduction in the spacing did not sig- 
nificantly change the singular values from the SVD computation Q. We considered the 
outputs to be a projection onto 4, 10 and 20 POD modes (corresponding to 4, 10 and 



20 mode output-projecti ons, as introduced in section 2.2.1 ). Using these modes, we use the 
expressions in equation ( 24|29 ) to obtain the matrices Ag/Bs, Cs defining the reduced-order 
model of the stable-subspace dynamics. The vorticity contours of the balancing and the 
adjoint modes, for a 10-mode output projected system, are plotted in Fig. [6j The adjoint 
modes provide a direction for projecting the linearized equations onto the subspace spanned 
by the balancing modes. Since these modes are quite different from the POD and the bal- 
ancing modes, the resulting models are also quite different from those obtained using the 
standard POD-Galerkin technique wherein an orthogonal projection is used. Since the mod- 
els obtained using balanced truncation are known to perform better than the POD-Galerkin 



models, as reported by Ilak & Rowley (2008), the better performance could be attributed 



to a better choice of projection using the adjoint modes. 

Since the reduced-order models of the stable-subspace dynamics are approximately bal- 



anced, the controllability and observability Gramians of the a^-dynamics of (24), given by 



expressions ([2]), are approximately equal and diagonal. Further, their diagonal values are 
approximately the same as the Hankel singular values obtained by the SVD The 
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Figure 8: Outputs (projection of the flow field onto POD modes) from a reduced-order 
model obtained using a 20-mode output projection. The first (top figure) and eleventh 

(bottom figure) outputs of the DNS ( , o) are compared with predictions of models 

with 4 ( , x), 10 ( , V), and 20 ( , ) modes. 



diagonal values of the Gramians and the singular values for different output projections are 
plotted in Fig. [7| for a 30-state reduced-order model. With increasing order of output pro- 
jection, the HSVs converge to the case with full-state output, and the number of converged 



HSVs is roughly equal to the order of output projection, as was observed by |Ilak fc Row- 



ley (2008). We see that the diagonal elements of both the Gramians are very close to the 
HSVs for the first 20 modes. For higher modes, the observability Gramians are inaccurate, 
which is due to a small inaccuracy of the adjoint formulation mentioned in section [STT] For 
controller design, we use models of order < 20, for which these Gramians are sufficiently 
accurate. 

Finally, to test the accuracy of the reduced-order models, we compare the impulse responses 



of system ( [44| ) (that is, restricted to the stable subspace) with that of the model (24), 
restricting au = 0. In particular, we compare the outputs of the two systems, which are 
the projection onto the POD modes; a representative case in Fig. [8] shows the results of 
4, 10 and 20 mode models of a system approximated using a 20-mode output projection 
(the outputs are projection onto the leading 10 POD modes). The first output, which is a 
projection onto the first POD mode, is well captured by all the models until t ^ 60, while 
the 20-mode model performs well for all time. Also shown is the eleventh output, which is 
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Figure 9: Schematic of the implementation of fuh-state feedback control in the nonlinear 
simulations. The entire velocity is first projected onto the unstable eigenvectors and the 
stable subspace POD modes to compute the reduced-order state a. The state is then 
multiplied by the gain computed based on the reduced-order model using LQR, to 
obtain the control input u. 
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Figure 10: Control gain obtained using LQR, and the initial condition is that obtained 
by an impulsive input to the system. Control is turned on at t = 0. Comparison of the 

outputs Hui and yss of a 12-mode reduced-order model ( , x) with the projection of data 

from the linearized simulation ( , o). 

well captured only by the 20-mode model. As we will see later, it is important to capture 
the higher-order outputs for design of observers. 



4.4 Full-state feedback control 



The resulting models can now be used along with standard linear control techniques to 
obtain stabihzing controllers. We u^e Linear Quadratic Regulator (LQR) to compute the 
gain K so that the eigenvalues of {A + BK) (where the matrices were defined in (24)) are 
in the left-half of the complex plane, and the input u = Ka minimizes the cost 



J[a, u] 



f 

Jo 



{a^Qa + u'Ru) dt^ 



(46) 



where Q and R are positive weights computed as follows. We choose Q such ^hat the first 
term in the integrand of (46) represents the energy, that is, we use Q — with C 
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Figure 11: Lift-coefficient Cl vs. time for full-state feedback control, with control turned 
on at different times in the base uncontrolled simulation. The base case with no control 

( ) has the unstable steady state as the initial condition, and transitions to periodic 

vortex shedding. The control is tested for different initial conditions, corresponding to t — 
170, 180, 210 of the base case, and stabilizes the steady state in all the cases ( ). 



defined in (25). The weight R is chosen to be a multiple of the identity c/, and typically c 
is chosen to be a large number ~ 10"^"^, to avoid excessively aggressive controllers. The 
control implementation steps are sketched in Fig.[9| first compute the reduced-order state a, 
using the expression (31), then the control input is given hy u — Ka. Here, we derive the 



gain K based on a 12-mode reduced-order model (with 2 unstable and 10 stable modes), 
using R — 10^, and include the same in the original linearized and nonlinear simulations. 
The output is approximated using a 4-mode output projection. The difference between the 
linear and nonlinear simulations is that, in the latter, the steady state field is subtracted 
from the state x, before projecting onto the modes to compute the reduced-order state a. 



Fig. [To] compares the model predictions with the projection of data from the simulations 
of the linearized system (38, 39), with a control input. The initial condition used is the 



flow field obtained from an impulsive input to the actuator. Both the states shown in the 
figure eventually decay to zero, which implies that the perturbations decay to zero, thus 
stabilizing the unstable steady state. More importantly, the model predicts the outputs 
accurately for the time horizon shown in the plots. 



We now use the same controller in the full nonlinear simulations and test the performance 
of the model for various perturbations of the steady state. A plot of the lift coefficient Cl 
vs. time t, with the control turned on at different times of the base simulation, is shown 
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Figure 12: Outputs of a system with full-state feedback control. The control gain is obtained 
using LQR, and the initial condition is that corresponding to t = 180 of the uncontrolled 
case plotted in Fig. [TTj Comparison of the outputs yui and ysi of a 12-mode (2 unstable 

and 10 stable modes) reduced-order model ( , x) with the projection of data from the 

full nonlinear simulation ( , o). 



in Fig. 11, The initial condition for the base case (no control) is the unstable steady 
state; eventually, small numerical errors excite the unstable modes and the flow transitions 
to periodic vortex shedding. In separate simulations, control is turned on at times t = 
170, 180, 210 corresponding to the base case. As the figure shows, the control is effective 
and is able to stabilize the steady state in each case, even when the flow exhibits strong 
vortex shedding. We remark that the latter two of these perturbations are large enough to 
be outside the range of validity of the linearized system, but the control is still effective, 
implying a large basin of attraction of the stabilized steady state. We also compare the 
output of the reduced-order model with the outputs of the nonlinear simulation; the plots 



are shown in Fig. [T2j The models perform well for the initial transients, but for longer 
times fail to capture the actual dynamics. This is not surprising as these perturbations 
are outside the range of validity of the linear models. For control purposes, it appears to 
be sufficient to capture the initial transients (approximately one period), during which the 
instability is suppressed to a great extent. We remark that one could possibly compute 
nonlinear models by projecting the full nonlinear equations onto the balancing modes, or 
enhance the model subspace by adding POD modes of vortex shedding and the shift modes 



as proposed by Noack et al (2005) to account for the nonlinear terms. 



Finally, we note that the reduced-order model (24) decouples the dynamics on the stable and 



unstable subspaces, and also, the dynamics on the unstable subspace can be computed only 
using the unstable eigen-bases ^tnd "i^u- Thus, we could derive a control gain K G R-^^^^, 
based only on the two-dimensional unstable part of the model, such that the eigenvalues 
of {Au — BuK) are in the left half complex plane. That is, we can obtain a stabilizing 
controller without modeling the stable subspace dynamics. We have performed simulations 
to test such a controller and found that it also is capable of suppressing the periodic vortex 
shedding and thus results in a large basin of attraction for the stabilized steady state. The 



choice of weight matrices Q and R in the LQR cost (46) needs to be different to obtain a 
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Figure 13: Schematic of the implementation of observer-based feedback control in the non- 
linear simulations. The control input u and the sensor measurements y are used as inputs 
to the observer, which reconstructs the reduced-order state a. This state is then multiplied 
by the gain to obtain the control input u. Both, the controller and observer gains K 
and L are computed based on the reduced-order model using LQR and LQG respectively. 

comparable performance. However, as shown in the next section, it is essential to model 
the stable subspace dynamics to design a practical controller based on an observer that 
reconstructs the entire flow field using a few sensor measurements. 



4.5 Observer-based feedback control 



The full-state feedback control of section |4.4| is not practical since it is not possible to 
measure the entire flow field. Here, we consider a more practical approach of measuring 
certain flow quantities at a small number of sensor locations. We assume that we can 
measure the velocities at the sensors shown in Fig. [ij in the near-wake of the plate. We 
remark that, even though these sensors are not experimentally realizable, they serve as a 
good testing ground for our models. 



4.5.1 Observer design 



Using the reduced-order models derived as outlined in section [431 we design observers that 
dynamically estimate the entire flow field. Since the models (24) have a different output 



(projection onto certain modes), we modify the output equation so that it appropriately 



represents the sensor measurements. We replace the output equation of (24) with 



def -F7 

a — Ga, 



(47) 



where M G and s is the number of sensor measurements. The matrix M is sparse 

and extracts the values of the output of ([24|) at the sensor locations; thus, each row of M is 
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filled with Os except for the entry corresponding to a sensor measurement, which is 1. With 



the output equation (47), we design an observer using a Linear Quadratic Gaussian (LQG) 
estimation. This method assumes that the errors in representing the state a and and the 
measurement y (due to the inaccuracies of the model) are stochastic Gaussian processes, and 
results in an estimate a of the state a that is optimal in the sense that it minimizes the mean 
of the squared error; refer to Skogestad & Postlethwaite (2005| for details. We now discuss 



briefly our procedure for modeling these noises; consider the reduced-order model (24), but 
with process noise w and sensor noise v which enter the dynamics as follows: 



a — Aa + Bu + w 
y — Ca + V. 



(48) 
(49) 



A key source of the process (state) noise w arises from model truncation, and second, from 
ignoring the nonlinear terms in the reduced-order model. The nonlinearity of the dynamics 
is important, for instance, when the model is used to suppress vortex shedding. A source 
of the sensor noise arises from two sources; first, the state x is approximated as a sum 



of a finite number of modes (17), and second, in the output projection step, the output 



is considered as a projection of the (approximated) state x onto a finite number of POD 



modes (22). Here, we approximate these two noises as Gaussian processes whose variances 
are 



^ — f ifimeas) ^^measi 



(50) 
(51) 



and E(-) gives the expected value. Here, /(•) is the operator obtained by projecting the non- 



linear Navier-Stokes equations (32) onto the balancing modes using the adjoint modes 
The state aj^eas is obtained by projecting the snapshots, obtained from a representative 
simulation of the full nonlinear system, onto the balancing modes. The representative sim- 
ulation we used here is the base case, with no control, shown in Fig. 11, which includes 



the transient evolution from the steady state to periodic vortex shedding, 
estimator is of the form: 

^ Aa + Bu + L{y - Co), 
y Ca, 



The resulting 

(52) 
(53) 



where a is the estimate of state a, y is the estimated output, and L is the observer gain. The 



estimator is then used along with the full-state feedback controller designed in section 4.4 



to determine the control input; a schematic is shown in Fig. 13 



Since the observability Gramian corresponding to the pair (A, C) is different from that 
for the pair (A, C), the model (24) with the output represented by (47), is not balanced. 
In principle, it is possible to construct reduced-order models of a system with the sensor 
measurements as the outputs, using the procedure of section |2.3[ Since the number of 
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of ([44|), and the energy captured by 4 ( , x), 10 ( , V), and 
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Figure 16: Lift-coefficient Cl vs. time t, for estimator-based feedback control, with control 
turned on at different times in the base uncontrolled simulation. The base case ( ) is the 



same as in Fig. [TT| and the control is tested for different initial conditions, corresponding 

to t = 170, 180, 210 of the base case ( ). In both the cases, the controller stabilizes the 

flow to a small neighborhood of the steady state. 



outputs will be typically small, the output projection step would not be required for such 



models. However, if this were done, the cost function (46), based on total kinetic energy in 
the perturbation, would not be captured well by the model and the model would not be as 
effective for control design. 

The models used in section [44] for full-state feedback were those of a system whose stable- 
subspace output was the velocity field projected onto the leading 4 POD modes. These 
4 POD modes capture only about 85% of the energy, but the resulting models were effective 
in suppressing vortex shedding. However, for observer design, this representation of the 
output is inadequate, as the energy content of the flow at the sensor locations is very small, 
while the POD modes capture the energetically dominant modes. Hence, a greater number 
of POD modes is required to accurately represent the velocity at the sensor locations. The 
temporal evolution of the energy content of the flow, obtained from an impulse response of 
the system restricted to evolve on the stable subspace, is plotted in Figure[T4| Also plotted is 
the energy content of the same flow, but projected onto the leading 4, 10 and 20 POD modes; 
thus, a 4-mode projection leads to noticeable errors, while both 10- and 20-mode projections 
accurately represent the energy. The velocity field at the sensor locations, reconstructed 
by 10 and 20 POD modes, plotted in Figure [15) shows that a 10-mode projection does 
not accurately represent the velocities at the sensor locations. Since 20 POD modes are 
sufficient to represent these velocities, we derive models using a 20-mode output projection, 
and use the same for observer design. 
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Figure 17: States of the system with observer-based control; the states recon- 
structed ( , x) by a 22-mode observer quickly converge to the actual states ( , o). 



The initial conditions used are those corresponding to t — 180, 210 (top and bottom) of the 



uncontrolled case shown in Fig. 16 



4.5.2 Observer-based control 



The models obtained using the modified output (47) are used to design dynamic observers 
based on the vertical {v-) velocity measurements at the sensor locations. A 22-mode 
reduced-order model, with 2 and 20 modes describing the dynamics on the unstable and 
stable subspaces respectively, is used to design a Kalman filter for producing an optimal es- 
timate of the velocity field based on Gaussian approximations of error terms (50, 51). This 
estimate is then used along with reduced-order model controller to determine the control 
input, as shown in Fig. [l3j The results of this observer-based controller (or compensator) 
are shown in Figs. 16, TTj The compensator again stabilizes the unstable operating point, 
and furthermore, the observer reconstructs the reduced-order model states accurately. Ini- 
tially, the observer has no information about the states (the initial condition is zero), but 
it quickly converges to and follows the actual states. There is a key difference from the 
full-state feedback control, that the compensator does not exactly stabihze the unstable 
operating point but converges to its small neighborhood. The reason is that the observer 
design is based on the velocities at the sensor locations projected onto the leading 20 POD 
modes, rather than the exact velocities at these locations. These small errors enter into the 
observer's dynamics in the same way that sensor noise enters, resulting in small errors in 
the state estimate. 
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5 Summary and discussion 



We presented an algorithm for developing reduced-order models of the input-output dy- 
namics of high-dimensional linear unstable systems, extending the approximate balanced 



truncation method developed by Rowley (2005) for stable systems. We assumed that the 



dimension of the unstable eigenspace is small and the corresponding global eigenmodes can 
be numerically computed. The modeling procedure treats the dynamics on the unstable 
subspace exactly and obtains a reduced-order model of the dynamics on the stable subspace. 

In a proof-of-concept study, the procedure was applied to control the 2-D low-Reynolds- 
number flow past a flat plate at a large angle of attack a, where the natural flow state 
is periodic vortex shedding. We first performed a continuation study at i?e = 100 and 
computed the branch of steady states with a varying from to 90°; we show that the 
flow undergoes a Hopf bifurcation from steady state to periodic shedding at a ^ 27°. 
We developed reduced-order models of the linearized dynamics dX a — 35° actuated by a 
localized body force close to the leading edge of the plate. The outputs were considered 
to be the entire flow field, projected onto the unstable eigenmodes and the leading POD 
modes of the impulse response simulation (restricted to the stable subspace). We developed 
stabilizing controllers based on the reduced-order models to stabilize the unstable steady 
state and showed that the models agreed well with the actual simulations. We also included 
the controllers in the full nonlinear simulations, and showed that they had a large-enough 
basin of attraction to even suppress the vortex shedding. For such large perturbations, 
however, the model agreement with the full simulation was good only for short times. A 
natural step towards improving these models would be to project the full nonlinear equations 
onto the balancing modes to obtain nonlinear models. Alternately, the balanced models, 
which accurately capture the transient dynamics, could be combined with the POD-based 



models using shift-modes of Noack et al. (2005) which accurately capture vortex shedding 



and some of the transient dynamics. An interesting future direction is development of 
algorithms to compute nonlinear balanced models, for instance based on the theoretical 



work of |Scherpen| ( |1993| ) 



Instead of computing nonlinear models, here we pursued a step towards more practical con- 
trollers by considering an observer-based control design, in which the outputs were modified 
to be just two near- wake velocity measurements. The nonlinear terms in the equations, 
which our models do not capture, were treated as process noise, and the error in modeling 
the outputs was treated as sensor noise. We designed a 22-mode reduced-order observer 
which reconstructed the flow field accurately, and along with the controllers, suppressed vor- 
tex shedding and stabilized the flow in a small neighborhood of the unstable steady state. 
We remark that the actuator and sensors considered here are not practically realizable, 
but the methodology presented here can be extended to a more practical actuation such as 
blowing and suction through the plate and measurements using surface pressure sensors. 
Furthermore, the choice of sensor locations in this study was ad hoc, and an interesting 
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problem is of finding the optimal sensor locations, for a given actuator. The controllers 
present here are designed to operate at a fixed set of parameters (such as Re and a), and it 
would be interesting to test their performance at off-design parameter values; the study on 



the performance of models of the linearized channel flow at off-design Re by |Ilak fc Rowley 
( [2008) shows promise in that direction. 

A motivation for the choice of our model problem was to develop tools towards manipulating 
wakes of micro-air vehicles. Recently, Taira & Colonius (2009&) performed a numerical 



study of flow past low-aspect-ratio plates, and a future direction we intend to undertake 
is to perform a detailed continuation study of the same flow to explore the existence and 
stabilization of high-lift unstable steady states in this 3-D flow. 
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A Balancing transformation for unstable systems 

Without loss of generality, a transformation T (and its inverse) that decouples the stable 
and unstable dynamics of ([T]) can be written as: 



A. 

where the columns of and Tg span the unstable and stable right eigenspaces of A, while 
the columns of Su and Sg span the unstable and stable left eigenspaces of A. Further, 



these matrices are scaled such that S^Tu = luu S^Tg = Ins- The transformation (54) 



decouples the dynamics of ([!]) as given in ^ with the various matrices defined as follows: 

— 5'*AT^, Bu — S^B, Cu — CTu, 
As = StATs, B, = S*B, Cs = CTs. (55) 

Using ( 54 ) in Q , the Gramians of the original system ([T]) are 

Wo^S*XSu + SlW',Ss, (56) 
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where, and are the Gramians corresponding to the system defined by (A^, Bg^Cg), 
while and are the Gramians corresponding to the system defined by {—Au, B^^ Cu)- 
Let ^ R^^^^^ be the transformation that balances the Gramians and W^^ while $5 G 
j^risxris ^j^g transformation that balances Vl/^^ and W^. Then, it can be verified that the 
transformation that balances the Gramians Wc and Wo is given by 

def 



(57) 



Thus, the balancing transformation consists of two parts and $s which respectively 
balance the dynamics on the unstable and stable subspaces of A. As per the technique 



of Zhou et al. (1999), a reduced-order model can be obtained by truncating the columns 
of $ that correspond to the relatively uncontrollable and unobservable states. As we will 
show now, the algorithm outlined in section |2.3| essentially computes the leading columns 
of $5 (and the corresponding rows of its inverse). We show that the controllability Gramian 



of the stable dynamics of dTl), which are defined by (15), is the same as the "stable" part of 
the Gramian defined in ( |56[ ). Note that using (54) and the definition (14), the projection 
operator Pg can be written as 



Using the definition ([2]), the controllability Gramian of (15) is 



(58) 



f 

Jo 

f 

Jo 



sAt( 



.B) 



^^)*^ dt 



^StATst 



S:B 5*5. e 



T:A*Sst 



T* dt using equation (58) 



^Ast 



BsB^ e 



Alt 



dt ]T: 



using equation ( 55 ) 



-^s^^ c-^s ' 



(59) 



which is the same as the stable part of Wc- Similarly, it can be shown that the observability 
Gramian of (15) is the same as the "stable" part of the observability Gramian Wq' 



roc 

/ e^^*^**(P:C*) (P:C*)*e(^^^*)** dt 
Jo 



StW'Ss. 



(60) 



Thus, balancing the Gramians and is identical to balancing the parts of the Grami- 
ans Wc and Wo of the original system ([T]) that are related to the dynamics on the stable 
subspace of A. 



B Derivation of the adjoint equations 



In this appendix, we derive the adjoint of the linearized semi-discrete equations ( 38|39 ). 



Let (C t/j) be the weighting functions corresponding to (7, /). Then, using the inner product 
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defined in equation (40), the weak form of ( 38|39 ) is: 



dx dt 



(61) 



+ / / i> ■ ECsdxdt = 0. 
Jo Ja 

Integrating by parts with respect to t and rearranging terms, 

£ lr{- (C^Cr'f^ + (C^Cr'C^E^i: + PC- {{C^C)-'C^AfUlo)fc) dxdt 



+ £ J^f. (^ECiC^C)-\) dx dt + (7, C) 



0. 



(62) 



For linearization about stable steady states, 7 ^ 0, as T ^ oc, and if the adjoint equations 
are integrated backwards in time, ({t = 0) ^ 0. So, the last term on the left hand side of 
equation (62) vanishes identically. If equation (62) is to hold for all values of 7 and /, we 
get the following adjoint equations hold: 



dt 



+ C'E'i, = -f3C' CC + C)Ml{ioY qa, 
ECi = 0, 



(63) 
(64) 



where ^ = {C^C)~^C, and qa — can be thought of as the weighting functions corre- 
sponding to the streamfunction s and the flux q respectively. Now, equations ( 63|64 ) have 
the same form as ( 38|39 ) except for the nonlinear term. Thus, the same time-integration 
scheme can be used for both, with the appropriate (linearized) nonlinear terms. 
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